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NOMENCLATURE 


Symbols 


Solar collector receiver surface area 

One half the solar collector aperature area 
Cantilevered beam width 

Solar CGOllector depth 

Specific heat 

Solar collecevor concentratiren ratio 


Young's modulus of elasticity, or the 
penalty function exponent used in DESOP 


Lower Golden Section fraction 
Upper Golden Section fraction 
Cantilevered beam height 
Equality constraint vector 


A user's flag in the DESOP program for 
initial and final user generated output 


Cantilevered beam length, or the solar 
collector length 


Solar collector mass flow rate 
Solar collector mass flow rate 
Number of design variables 
Objective function 

Penalized objective function 
Cantilevered beam load 


Solar collector heat available 
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cae = 
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Thetal = 
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Solar collector heat gain 
Yearly solar collector heat available 
Heat flux 


Penalty parameter used in DESOP or the 
Solar collector receiver radius 


Solar collector receiver radius 

Solar collector aperature cover temperature 
Solar collector coolant exit temperature 
Solar collector acceptance half angle 

Solar collector truncation angle 


Solar collector distance between the 
reflector and a point tangent to the receiver 


Cantilevered beam volume 
Solar collector reflector coordinate 
Design variable vector 


Solar collector reflector coordinate 


Greek Letter Symbols 


Qa _ 


One-dimensional search move parameter 
Cantilevered beam deflection 


Solar collector geometry angle measured 
from the collector centerline (See Fig. 6) 


Solar collector acceptance half angle 
Solar collector truncation angle 
Cantilevered beam shear stress 


Cantilevered beam bending stress 
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A. BACKGROUND 

Most engineering design problems contain several con- 
tinuous variables and as such have an infinite number of 
solutions. The purpose of optimization is to find the best 
possible solution among the many potential solutions for 
a given problem in terms of some effectiveness or performance 
criteria. There are several methods of optimization. The 
methods may be classified as follows: 

Analytical methods which use the classical techniques 
of differential calculus and the calculus of variations. 

Numerical methods which use past information to generate 
better solutions to the optimization problem by means of 
1terative procedures. Numerical methods can be used to 
solve problems that cannot be solved analytically. 

Saopmical methods which use the preparation of a plot 
of the parameter to be optimized as a function of one or 
more variables. This method although simple and easy to 
use becomes unmanageable when there are three or more 
design variables. 

Experimental methods which use direct experimentation 
of the actual process, the results of one experiment being 


used to decide on where to perform the next experiment. 


iz 





Case study methods which evaluate the results from a 
number of representative cases, and choose the "best" 
solution. The "best" solution is thus not likely to be 
the optimum solution. 

Of the optimization methods, the numerical method lends 
itself to computerized solution. As design is an inter- 
active process between the designer and the problem, and 
the desktop computer lends itself towards dedicated inter- 
active use, the development of a numerical optimization 
program for use on a desktop computer in an interactive 
mode, is the objective of this thesis. 

To date a great amount of effort has been spent de- 
veloping reliable and efficient optimization programs for 
mainframe computers. These programs are fairly large and 
complex, requiring a substantial amount of core space 
during execution. The size and complexity of the programs 
has been the result of an attempt to minimize the amount 
of computer time required to perform an optimization and 
thus the cost to the user. With the advent and availability 
of desktop computers, there has been a sharp reduction in 
the cost of computer time to the computer user. While the 
desktop computer has far less core space than a mainframe 
computer, once the time factor is removed from the numerical 
optimization process it is possible to put a small but 
reliable numerical optimization program on a desktop 


computer. A design problem concerning the optimal geometry 


iS 





for a nonimaging parabolic trough solar collector was 
developed to demonstrate the engineering application of the 


numerical optimization program developed for this thesis. 


Pee oCOPE 

The numerical optimization of a given function may be 
accomplished using many varied and different algorithms. 
Some of the more popular methods used are: random search, 
linear programming, feasible directions, Golden Section, 
Newton's method and sequential unconstrained minimization. 
A particular optimization program will use one or more of 
these methods to efficiently and reliably arrive at the 
best solution to a particular problem. 

For this thesis two computer programs were developed 
to perform numerical optimization on a desktop computer. 
The first program was developed as an executive program to 
control the optimization process. The executive program 
is named OPCON which stands for OPtimizer CONtrol program. 
OPCON provides the interactive link between the user and 
the program which performs the numerical optimization. 
OPCON allows the user to input data, attach a specified 
analysis subprogram to the numerical optimization program 
and control execution of the numerical optimization program. 
The second program developed was the numerical optimization 
program, DESOP. DESOP stands for DEsktop Sequential uncon- 


strained minimization technique Optimization Program. 
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DESOP performs the numerical optimization of the user 
supplied problem using the seguential unconstrained mini- 
mization technique with an external penalty function. The 
unconstrained subproblem is solved using the Fletcher- 
Reeves method of conjugate directions, Golden Section 
search, and polynomial interpolation. 

The third computer program, NISCO, was developed to 
model a nonimaging concentrating compound parabolic trough 
solar collector using thermophysical, geophysical, optical 
and economic analysis to compute a life-cycle cost fora 
design with a stated energy capacity. NISCO stands for 
NonImaging concentrating compound parabolic trough Solar 


Beriector. 


See CbJECTIVE 

The objective of this thesis was to develop a system 
of interactive programs for the Hewlett-Packard 9845A 
desktop computer which perform numerical optimization, and 
to demonstrate the capability on the design of a nonimaging 
concentrating compound parabolic trough solar collector. 
Three programs were developed to meet the objective: an 
executive program, a numerical optimization program and 
a sOlar collector analysis program. 

The purposes of the executive program OPCON are: 

Memo provide a primary point of contact for the com-— 
puter user from which to effect a numerical optimization 


on any number of user prepared analysis subprograms. 
1b. 





2. TO provide a standardized, formatted input for the 
design variables, side constraints and optimizer control 
parameters, which is recognizable by all the programs in 
the numerical optimization package of programs. 

3. To control the operation of the different optimiza- 
tion and design analysis programs within the system through’ 
a process of program overlays which maximizes the computer 
Space available for the design analysis Progr al. 

4. To develop a program which 1s portable to different 
computers uSing an advanced BASIC language. 

The purposes of the numerical optimization program 
DESOP are to develop an optimization program that: 

1. Is relialbe in reaching a deSign optimum, irrespective 
@Seeene Starting point. 

2. Will arrive near the design optimum using default 
values for the optimizer control variables. 

3. Will allow the user to monitor the optimization 
process and to change the optimizer control variables to 
more efficiently and/or more accurately reach the design 
optimum. 

4. Is portable to different computers uSing an advanced 
BASIC language. 

The purposes of the solar collector program NISCO are 


iO: 
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1. Model a nonimaging concentrating compound parabolic 
trough solar collector using a system of thermophysical, 
geophysical, optical and economic equations. 

2. Arrive at an optimum design for a solar collector 
given a stated average daily heat gain and a life-cycle 


period. 
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rt UME tex OPT i Za ON 


A. THE CONCEPT OF NUMERICAL OPTIMIZATION 

Consider the problem of designing the cantilevered beam 
shown in Figure 1. The design task may be broken down into 
three major parts. First, the objective of the design must 
be determined, which in this case is to minimize the weight 
of the beam required to support the concentrated tip load P. 
Second, any physical constraints that may effect the design 
must be determined. Thirdly, any limits which exist on the 
design variables must be determined. The design problem 
may then be reduced to a system of equations as follows: 

Minimize the volume (V) 

V = B-H-L 
Subject to: 


Bending stress (o)) 


Geb 
Oo, = “pep > 20000 psi 


Shear stress (v) 


oe 


V = 57B-H - 10000 psi 





Deflection under load (9) 


3 
aia; £ a aneh 
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P=10000 


oe 200 ee 


Pagure 1. 


Cantilevered Beam Design Problem 


19 





With geometric constraints such that: 


OMS. 7.0.5 
Oe soe 2.0. 0 
Byes) 10.0 


At this time the following definitions are introduced: 

Objective Function - The parameter which is to be 
minimized or maximized during optimization. The objective 
function always occurs on the left side of the equation 
unless it is also used as a design variable. An objective 
function may be either linear or nonlinear, implicit or 
explicit, but must be a function of the design variables. 

Design Variable - Any parameter which the optimiza- 
tion process is allowed to change in order to improve the 
design. Design variables appear only on the right hand 
Side of equations in the analysis program. 

Inequality Constraint ~- Any parameter which must 
not exceed specified bounds for the design to be acceptable. 
Constraint functions always appear on the left side of 
equations. A constraint may be linear, nonlinear, implicit 
or explicit, but must be a function of one or more design 
variables. 

Equality Constraint ~ Any parameter which must equal 
a specified value for the design to be acceptable. The 
same rules apply to equality constraints as inequality 


constraints. 


20 





Side Constraint - Any upper or lower bound placed 
upon a deSign variable. Side constraints are usually not 
included in the system of equations that comprise an 
analysis program. Instead they are usually included as 
part of the data input to the optimization program. 

Analysis Code - The system of equations utilizing 
the design variables which are used to calculate the objective 
function and the constraints of a particular design problem. 

The general optimization problem may thus be stated 
mathematically as: 

Find the set of deSign variables x. where 1 = 1,2,...,n 
which will: 

Minimize the objective function (Obj) 

Obj = £(X) 

Subj ect to: 

Inequality constraints (G) 


Come ores 0 ee eee 


Equality constraints (H) 
H; (X) = 20 ete eee, 


Side constraints 


xt <x. < xt Se cs ni 
ale Ey ale 


Returning to the cantilevered beam problem, it may be 
stated in the standard format as follows: 
feo <(l)o= s) X(2)— oH, and Obj >= Vol = B-H-L 


Then minimize Obj = Vol 


Zl 
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Memalezea objective function at the initial design may be 
found by taking the first partial derivative of Opj with 
respect to the deSign variables. The gradient of the 
penalized objective function defines the direction of 
steepest ascent. In the case of the cantilevered beam, it 
1S desired to minimize the objective function; therefore, 
the greatest improvement in design may be achieved by moving 
in the negative gradient, or steepest descent direction. 
From the initial deSign point a search is performed in the 
steepest descent direction for the minimum value of Opj in 
that direction. At the new minimum, the gradient of the 
penalized objective function is again determined anda 
search is performed in a conjugate direction until a second 
minimum is found. Successive iterations are performed 
until the gradient is found to be zero or each successive 
iteration produces a sufficiently small change in Op) such 
that for all practical purposes the minimum has been found. 
At this time the penalty function is increased. If the 
design is in a region where there are no constraints 
violated an increase in the penalty function will not 
change the value of Opj. If on the other hand the design 
is in an infeasible region where there are one or more 
constraints violated, Opj will be increased, and the 
search for a new minimum will commence. If the minimum 


of the objective function exists in the infeasible region, 
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the minimum value for the objective function in the feasible 
region will be approached from the infeasible region as the 
penalty function is increased. The design improvement 
process will terminate when a zero gradient is found or 
Successive iterations produce a sufficiently small change 
in the value of Opj and an increase in the penalty function 
causes no change in Opj. 

The minimum thus reached by the optimization process is 
a minimum with respect to the penalized objective surface 
immediately surrounding the final design point. The opti- 
mization process cannot distinguish between local and global 
minimum points. It is thus good engineering practice to 
run several optimizations for a particular design problem 
from several different initial design points. If optimiza- 
tions performed from different initial design points converge 
on the same minimum point, that point 1s probably a global 
minimum. If on the other hand two or more minimums are 
found, there may be local minimums located in the design 
Space being considered and care must be taken to find the 


global minimum. 


B. THE DESOP NUMERICAL OPTIMIZATION PROGRAM 

The DEsktop Sequential unconstrained minimization 
technique Optimization Program was developed using the 
basic optimization approach outlined in Section II.A. 


A copy of the program is included in Appendix E. The major 
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program structure is shown in Figure 2. The following 
discussion will refer to Figure 2 and describe the major 
features of the program. 

l. BaSic Program Execution 

The DESOP program begins execution when it is 

loaded, linked to the user's analysis subprogram, and the 
program is instructed to run by the OPCON program. The 
above actions are automatically performed by the OPCON 
program. DESOP is loaded into the computer by an overlay 
process. Therefore no variables can be directly transferred 
between the DESOP and OPCON programs. DESOP begins execu- 
tion by reading the optimizer control variables and the 
design variables that were input using the OPCON program 
and saved to a mass storage device. The program then sets 
Icalc equal to one and evaluates the objective function 
and constraints at the initial deSign point. Icalc 1s a 
flag provided the user to key user specified output on the 
initial and final design analysis. DESOP will provide the 
user with a hard copy output of the initial design variables, 
the value of the constraints, the objective function and 
the penalized objective function. The user then has the 
option of continuing with the DESOP program to optimize 
his analysis subprogram or to return to the OPCON program 


to change one or more of the input parameters. 
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Figure 2. DESOP Flow Diagram 





Proceeding with the optimization, there are two 
major loops in the optimization program. The outer loop 
increases the penalty function when the inner loop's con- 
vergence criteria have been met. A convergence test is 
then performed by the outer loop. If the convergence 
criteria is met, the optimization process is considered 
Finished. If the convergence criteria for the outer loop 
1s not met, program execution is returned to the inner 
loop. The inner loop performs successive iterations 
searching for the minimum of the penalized objective 
function with no increase in the penalty function taking 
place. When the inner loop's convergence criteria have 
been met program execution is transferred to the outer loop. 

Execution of the program while in the inner loop 
proceeds as follows: First, the gradient of the penalized 
objective function is calculated by subroutine GRAD. The 
program then computes a search direction using either a 
steepest descent method or the method of conjugate directions 
developed by Fletcher and Reeves [Ref. 1]. Once a search 
direction is established the optimizer attempts to locate 
the minimum value of the penalized objective function in 
the search direction. This process if referred to as the 
one-dimensional search and is illustrated in Figure 3. 

The efficiency and accuracy to which the one-dimensional 


search for the minimum of the penalized objective function 
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1s accomplished, is the key to successful sequential un- 
constrained minimization technique numerical optimization. 
The one-dimensional problem may be expressed in terms of 
the penalized objective function, Opj, and the amount of 
movement, a, in the search direction, S. First the slope 
of Opj with respect to a is calculated. An initial "guess" 
of how far to move 1S made uSing subroutine ALPGES. The a 
which corresponds to the minimum value of Opj in the one- 
dimensional search is then calculated using subroutine 
ALPBND and subroutine QUEBIC. The minimum value of Op3j 
thus found is then compared to the previous value of Op} 
for convergence uSing subroutine CONVRG. If convergence 
is not met, execution returns to the start of the inner 
loop. If convergence is met, execution returns to the 
outer loop. 

When the convergence criteria have been met for 
both the inner and outer loops, the program proceeds to 
set Icalc to three as a flag for user generated output for 
the final design. DESOP then provides the user with a 
hard copy output of the final design variables, objective 
function value, penalized objective value, constraint 
values, the number of inner loop iterations, the number of 
times the analysis subprogram was called and the final value 
Of the penalty function. The OPCON program is then overlayed 
Over the DESOP program and program execution is returned to 


the OPCON program. 
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2. Finding the Search Direction 
The first step in finding the search direction, S, 
is to determine the slope of Opj at the present design 


point. The forward finite difference method is used.where: 


oF _ F(X, - AX, ) - F(X, ) ae 
BK AX. a 
1 1 
tego e i 2, aie IGLay 
As JE/OX, gives the direction of positive slope, the search 
direction is the negative of OF/SX.. The first search is 


performed using the steepest descent as found above using 


the following relation: 


where alpha is the distance moved in the S direction. When 
a minimum is obtained along the direction of steepest 
descent, a new Fletcher-Reeves conjugate search direction 
[Ref. 1] is calculated at the new Design point using the 


following relations: 


Ss = ax + BS; 
ar 
Ndv yp! } 
BE [SS] 
: Ke 
— 1=1 a 
Ndv oF ,2 
E x. 
i=1 


where the prime denotes values for the present iteration 


and the non-prime variables indicate values for the previous 
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iteration. A one-dimensional search is then performed in 
the new search direction. Searches are continued using 
the conjugate direction method for Ndv + 1 iterations, 
where Ndv is the number of design variables. The search 
process is then restarted uSing the steepest descent method. 
The reason for incorporating the conjugate search method 
1s that the steepest descent method when traversing a design 
surface with a curved valley will tend to zigzag from one 
side of the design surface valley to the other making very 
little progress as is illustrated in Figure 4. The conju- 
gate direction method is much more efficient in traversing 
such a design surface. However, as the conjugate direction 
method is additive upon previous searches, it has a tendency 
to decrease in effectiveness with each successive search 
owing to the accumulation of numerical "noise." For that 
reason the search process is restarted with the steepest 
descent method every Ndv + 1 iterations, or when the con- 
jJugate direction predicts a positive slope. The search 
direction is normalized to avoid inaccuracies caused by 
numerical ill-conditioning. 
fee estlimacing am Initial Value for Alpha 

The initial estimate for alpha is made in the 
following manner: First, the slope of Opj in the search 
direction is calculated as the sum of each of the products 
of the gradients times the search direction. Then the 


Slope of Opj in the search direction is divided by the value 
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of Opj. This value is then multiplied by an improvement 
percentage in Opj. This first estimate is then applied to 
a series of conditional tests to determine the validity of 
the estimate with respect to the slope of Opj and the magni- 
tude of the design variables. Lastly, the estimate for 
alpha is checked to see if it violates any side constraints. 
If it does, the value of the estimate for alpha is reduced 
until the side constraints are no longer violated. 
4. Calculating Alpha 

Micemcaloularien @L alpha is the most critical al- 
gorithm in the DESOP program in providing reliable optimizer 
operation. The ability to accurately and efficiently find 
the minimum of the penalized objective function in the one- 
dimensional search affects directly the operation of the 
optimizer. Figure 5 illustrates the Zigzag phenomenon 
Which occurs if alpha is not accurately found. The Zigzag 
phenomenon is caused by the fact that the optimizer in per- 
forming the forward finite difference for calculating the 
search direction perturbs the design vector a very small 


amount. As such the optimizer can only "see" the design 
surface that is immediately adjacent to the design point. 
Therefore, if the minimum is not found in the one-dimensional 
search, the optimizer will converge very slowly on the 
minimum. 

There are two major sections to the ALPBND subroutine. 


The first section attempts to find the minimum value of Opj 
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Figure 5. The Zigzag Phenomenon 





uSing an expanding search technique. The first move is 
the amount predicted by the ALPGES subroutine. If the 
minimum is not bracketed by the first move, additional 
moves are made. Each move is larger than the last. The 
Size of the move is increased each time by an amount equal 
to the size of the last move divided by the lower Golden 


Section fraction, where the Golden Section fractions are: 


The lower Golden Section fraction, F 1s used so that the 


na 
interval will be consistent with the Golden Section search 
in the second section of the ALPBND subroutine. The expand- 
ing search is continued until the minimum value of the ob- 
jective function has been bracketed. 

Once the minimum is bracketed, a Golden Section 
search 1s performed to reduce the bracketing interval on 
the minimum by an amount such that when the two end points 
of the interval are taken with two points internal to the 
interval and a cubic is passed through the four points, the 
cubic will accurately predict the minimum of the penalized 
objective function. Himmelblau in [Ref. 2] states that the 
Golden Section search method of reducing the interval 


around the minimum of Opj is the most effective of the 


reducing techniques studied. Golden Section search is based 
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on the splitting of a line into two segments Known in 
ancient times as the "Golden Section." The ratio of the 
whole line ‘ea the larger segment is the same as the ratio 
of the larger segment to the smaller segment. The two 
Golden Section fractions are employed to split the interval 
bracketing the minimum as shown in Figure 3. Once the 
interval has been split, the two values of Opj corresponding 
to the internal points are compared to find the larger of 
the two. The internal point with the larger value of Opj 
will become the new end point for the interval, the remain- 
ing interior point will by the fact that it was determined 
by a Golden Section fraction, be equal to the point deter- 
mined by the other Golden Section fraction. Thus, only one 
new point must be calculated to continue the Golden Section 
Search. The search is continued in this manner until the 
vertical separation of the two end points with respect to 
the interior points is less than one percent. The four 
values of the penalized objective function corresponding to 
the four Golden Section search points are then sent to a 
BapLtc interpolator. The cubic interpolator will return a 
value for alpha to predict the minimum of the penalized 
objective function, and the minimum of the cubic function 
that the interpolator has created. The subroutine ALPBND 
will then test the predicted minimum with the minimum found 


at the predicted alpha. If there is less than a tenth of 
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one percent difference between the two values of the objec- 
tive function, the point predicted by the cubic interpolator 
will be accepted as the minimum and program execution will 
return to the main program. If the predicted minimum is not 
sufficiently close to the minimum at the predicted alpha, 
another Golden Section search will be performed to reduce 
the interval and better localize the minimum. The four 
points from the reduced interval will then be sent to the 
cubic interpolation subroutine. This process will continue 
until either the test for the minimum is positive or the 
interval has been reduced to less than 1LE-12. Program exe- 
cution will then return to the main program. 
pee subroutine QUEBIC 

Subroutine QUEBIC is used to estimate the alpha at 
which Opj iS a minimum based on four point cubic interpola- 
tion. If the function more closely resembles a quadratic 
than a cubic, a three point quadratic interpolation is 
performed using the three points which bracket the minimum. 
If the predicted minimum is outside the interval spanned 
by the two end points again a quadratic interpolation is 
performed. If the minimum still lies outside the two end 
points, the analysis returns to subroutine ALPBND, the 
inverval bracketing the minimum is reduced, and program 


execution returns to QUEBIC. 
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Seo ve ncencen or seae Penalized Objective Function 


The penalized objective function is tested for con- 
vergence at the end of each inner loop and again at the end 
of the outer loop in the main program. Convergence is 
tested by calling subroutine CONVRG. There are two criteria 
used for testing for convergence. The first tests the rela- 
tive difference of the value of Opj from the present iteration 
with the value of Opj from the last iteration. The second 
method tests the absolute difference of the two values. The 
second method is employed for cases when the value of the 
penalized objective function approaches zero. When conver- 
gence has been met on two succesSive iterations, the penalty 
function 1S increased by an amount specified by the user in 
the executive OPCON program. The penalized objective function 
is again tested for convergence. If convergence is still 
met, the optimizer considers the present value of the pena- 
lized objective function to be a minimum, noting again that 
numerical optimization programs cannot differentiate between 
local and global minimums. 

7. The Penalty Function 

The purpose of the penalty function is to increase 
the value of the objective function when the design is in 
an infeasible region. The infeasible region is that region 
where one or more design constraints are violated. When a 


constraint is violated, the value of the particular constraint, 
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aa is positive. The objective function is then penalized 


as follows: 


Opj = Obj + R-G,* 
where: 
Rea isd plueatien —COnseant 
E - an exponent constant 
This type of penalty function, one where the penalty is 
applied after the design leaves the feasible region, is 
known as an exterior penalty function. The exterior type 
of penalty function was chosen over other types, such as 
the interior or extended interior penalty function. If a 
function 1s discontinuous within the design space being 
studied, numerical difficulties may be encountered which 


make performing an optimization of the design difficult. 
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Pere COLLAR COLLECTOR OPTIMIZATION 


Peele CONCEPT OF NONIMAGING SOLAR COLLECTORS 

At the present time there are numerous schemes for 
the collection of solar energy and its conversion to a more 
useful form of energy. In the field of solar collectors 
there are three broad categories: flat plate collectors, 
focusing collectors and nonimaging collectors. The advan- 
tages and disadvantages of each type are shown in Table I. 
Welford and Winston in [Ref. 3] report that "in the mid- 
1960's, it was realized in at least three different labora- 
tories that light could be collected and concentrated for 
many purposes, including solar energy, more efficiently by 
nonimaging optical systems than by conventional image 
forming systems. The methodology of designing optimized 
nonimaging systems differs radically from conventional 
optical design. The new collectors approach very closely 
the maximum theoretical concentration; and for two-dimensional 
geometry, which is important for solar energy collection, 


il 


this limit is actually reached." Figure 6 shows the basic 


geometry for the nonimaging compound parabolic concentrating 


twelford, W. T. and Winston, R., The Optics of Nonimaging 
Concentrators, Light and Solar Energy, p. ix, Academic Press, 
mo 7 3 . 
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megure 6. 





Nonimaging Concentrating Solar Collector 
Geometry 
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mien collector depth, Cd, 1s used in the economic analysis 


Semcollector cost. The truncation angle, 6 is the maximum 


te! 
86 used in determining the collector geometry for the trun- 
cated collector. The collector may be significantly trun- 
cated before any appreciable change in the concentration 
ratio is affected. This allows a savings in manufacturing 
costs with little degradation in collector performance. 

The nonimaging concentrating solar collector will 
accept and deliver to the absorber all incident radiation 
that falls on the collector aperature and that is within 
the maximum acceptance half angles, Cs. That 1s, there is 
eaearcc of sky, 20., Paoniwrken  aMieeradratien,. Glrect, difiuse 
and reflected, is delivered to the absorber. It is this 
fact which makes the nonimaging concentrating collector 
attractive for solar energy use. Figures 7, 8, and 9 show 
ray paths for a nonimaging, truncated collector for the 
following cases: In Figure 7 the solar altitude is within 
the acceptance half angles. In Figure 8 the solar altitude 
1S equal to the acceptance half angle. In Figure 9 the 
solar angle is less than the acceptance half angle. [In 
Figure 7, all the radiationthat enters the collector is 
delivered to the absorber tube and is somewhat scattered 
over the absorber surface. In Figure 8 again all the 
radiation that enters the collector is delivered to the 
absorber tube, but is now tangent to the tube and 1s con- 


centrated on the front edge of the absorber tube. In 
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ieure 7. Ray Path Drawing for Solar Altitude Within the 
Maximum Acceptance Half Angle 
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Figure 8. Ray Path Drawing for Solar Altitude Equal to 
the Maximum Acceptance Half Angle 
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Figure 9 all the reflected radiation leaves the collector, 
thus the acceptance half angle provides a very sharp cutoff 


angle for accepting incident radiation. 


Bee SOLAR COLLECTOR DESIGN PROGRAM 

@& Nonimaging compound parabolic trough Solar COllector 
hereafter referred to as NISCO was chosen to model in the 
SeysiS program. There are seven major heat flow paths 
considered in the program. The major heat flow paths are 
shown in Figure 10 and are outlined below: 

q, 7 The sum of the direct, diffuse and ground reflected 
radiation that is incident on the collector cover. 

qd, ~ The sum of the direct, diffuse and ground reflected 
radiation that is reflected by the collector cover. 

a3 7 The sum of the direct, diffuse and ground reflected 
radiation that is absorbed by the collector cover. 

ah. The sum of the direct, diffuse, and ground reflected 
radiation that 1s transmitted by the cover and delivered 
either directly or indirectly to the absorber and is absorbed 
by the absorber. 

qs ~ The sum of the radiation reflected by the absorber 
that is absorbed by the cover and the thermal radiative and 
convective exchanges between the absorber and the cover. 

de - The sum of the thermal radiative and the convective 
exchanges between the cover and the environment. 


qz- The energy delivered to the collector coolant. 
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Penmaes O. Basie Solar Collector Heat Paths 
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To determine the energy balance on the collector, a 
procedure outlined by Kreith and Kreider, [Ref. 6] is followed. 
All important heat fluxes are first calculated from basic 
heat-transfer principles. The fluxes are then combined in 
heat-balance equations for the receiver, the aperature cover, 
and the coolant fluid. Since the various flux terms are 
nonlinear in temperature a Simultaneous iterative solution 
1s used to solve the equations. Note that the first page of 
Appendix Gis a cross reference list of the major equations 
used in the subprogram NISCO and the sources used to obtain 
the equations. 

The following discussion will detail the procedure used 
in the NISCO subprogram to calculate the heat gain for a 
particular solar collector design and the resultant life- 
cycle fuel savings. Variable names and program line numbers 
correspond to those used in the NISCO subprogram. The design 
variables chosen for the NISCO subprogram were: (1) Thetal, 
(O.), the maximum acceptance half angle, (2) Thetat, (8,), 
the truncation angle, (3) R, (r), the receiver radius, (4) L, 
ime collector length, and (5) Mfr, (m), the coolant mass flow 
rate. The objective function of the NISCO subprogram is the 
life-cycle fuel savings of the collector. Life-cycle fuel 
Savings are calculated as the cost savings realized by the 
collector over purchasing natural gas per unit quantity of 
10° Btu's expressed in present worth using the present worth 


analysis described by Newnan in [Ref. 13]. The constraints 


50 





placed upon the design by the NISCO subprogram were: the 
maximum and minimum values allowed for the truncation angle 
due to collector geometry as specified in Section III.A., 
the mass flow rate of the coolant must be positive, the 
receiver radius must be positive, a minimum average daily 
heat gain, and a maximum allowable coolant temperature. 
The side constraints placed upon the design to ensure that 
the final design was a reasonable design were: the maximum 
acceptance half angle must be greater than three degrees but 
less than 85 degrees, the truncation angle must be greater 
than 185 degrees but less than 260 degrees, the receiver 
radius must be greater than one tenth inch but less than two 
inches, the collector length must be greater than five feet 
but less than 100 feet, and the mass coolant flow rate must 
be greater than five lbm/hr and less than 1000 lbm/hr. 

The NISCO analysis subprogram proceeds as follows. 
A minimum average daily heat gain is specified and entered 
as Ql in line 1655. The subprogram is then SAVED and the 
OPCON/DESOP optimization program run. When DESOP calls 
the NISCO subprogram, it will pass the design variable 
vector X to the NISCO subprogram. The NISCO subprogram in 
lines 205 to 225 sets the design variable vector X equal 
to the design variables used in the NISCO subprogram. 
NISCO then proceeds to read in the constants and data used 


in the solar collector design. The design specifications 


a. 





are summarized in Appendix B. Collector geometry calcula- 
tions are then performed in lines 820 to 855 to ascertain 
the optical properties of the collector as developed in 
[Ref. 3]. An initial receiver and aperature cover tempera- 
ture iS assumed in lines 865 and 870. Monthly calculations 
are then performed to calculate the collector tilt angle, 
the minimum accepted solar altitude, the ground angle factor 
and sky heat loss constants as specified in [Refs. 5 and 6]. 
An hourly calculation is then performed to calculate the 
angle that the sun makes with the collector aperature cover, 
the amount of radiation incident on the collector aperature 
cover and the collector cover transmissivity and absorptance 
as specified in [Refs. 5 and 6]. The iterative portion of 
the analysis then proceeds as follows: First, the heat 
transfer coefficients for the collector are calculated as 
prescribed in [Refs. 4, 6 and ll], lines 1060 to 1145. Next 
a heat balance is performed on the cover as prescribed in 
[Ref. 6] and a new collector aperature cover, Tap, is cal- 
culated in lines 1155 to 1210. A heat balance is then 
performed on the receiver as specified in [Ref. 6] and the 
energy passed to the coolant, Qu, is calculated in lines 
1225 to 1305. Finally a heat balance is performed on the 
coolant as prescribed in [Ref. 6] and the coolant exit 
temperature is calculated in lines 1320 to 1365. Knowing 


the coolant exit temperature, Tc2, an average receiver 


ae 





temperature may be calculated and compared to the initial 
assumed receiver temperature, lines 1375 to 1405. If the 
new receiver temperature is within a tenth of one percent 

of the old receiver temperature, the analysis program 
proceeds to calculate the pumping power required for the 
collector. If the new receiver temperature 1S not eens 
the convergence specified, the iterative heat balance 
process is repeated, setting the aperature cover temperature 
and the receiver temperature to the new value calculated. 
The pressure drop through the collector and the power 
required to pump the coolant through the collector are 
calculated in lines 1430 to 1480, as specified in [Refs. 6 
Poe?) . The energy required to pump the coolant through 
the collector is then subtracted from the energy gained by 
the collector to calculate the available collector energy, 
Qa, line 1490. Qa is then summed for each hour that the 
analysis is performed. As the analysis is performed for 

one day each month, the summation of the available collector 
energy is then multiplied by thirty to obtain a yearly heat 
gain, Oy, line 1590. The cost of an equivalent amount of 
natural gas is then calculated. Using present worth analysis 
as described in [Ref. 13], a life-cycle savings is then 
calculated and the initial manufacturing cost of the 
collector is calculated and subtracted from the life-cycle 


Savings. This result is then divided by the life-cycle 


oe 





useful energy gain by the collector and multiplied by 10° 

to obtain the life-cycle fuel savings used as the objective 
function for the optimization program. The constraints on 
the design are then calculated in lines 1660 to 1685. Lines 
Moto 1815 contain printout specifications for the first 
and the last time that the NISCO subprogram is called by 

the optimization program. The values of the objective 
function and the constraints are then passed along with 


@eegram Control back to the optimization program. 
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m6 CURESULATS 


A major portion of the thesis work was devoted to the 
development of a reliable optimization program that would 
optimize a wide variety of problems. Following the accom- 
plishment of this goal, a subprogram was developed to 
model a nonimaging concentrating compound parabolic trough 


solar collector. 


A. RESULTS OF THE DESOP PROGRAM DEVELOPMENT 

In developing the DESOP numerical optimization program, 
four standard numerical optimization test problems were 
used. The four test problems provided a wide variety of 
different numerical problems with which the optimization 
program had to deal. The four test problems are listed in 
Appendix F. A major goal of this thesis was to optimize 
all four problems using default optimizer control variables. 
The results of the DESOP program optimizing the four test 
problems are given in Appendix A. The DESOP program was 
able to make significant design improvements in all four 
test problems uSing default optimizer control parameters. 
When the optimizer control parameters were adjusted for 
each individual problem, the DESOP program's performance 
was improved in all four cases. Further experimenting 
with the optimizer control parameters could lead to an 


even better performance of the DESOP program. 
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Pee RESULTS OF THE NISCO SUBPROGRAM 

A listing of the NISCO subprogram is included as 
Appendix G. There are numerous comment statements in the 
subprogram. The reader 1s encouraged to look closely at 
the subprogram. The results of the NISCO subprogram may 
be found in Appendix B. The first section of Appendix B 
details the design specifications that were chosen for the 
solar collector model. Three different daily heat capa- 
cities were specified 10000, 30000 and 50000 Btu, and the 
NISCO subprogram was used to find the optimum design for 
each. The second section of Appendix B gives the initial 
design and final designs for the three solar collector 
capacities. In each case the DESOP program was able to 
Significantly improve the design. The instantaneous 
efficiencies for the final designs are within a few 
percent of the instantaneous efficiency reported by Kreith 
and Kreider [Ref. 6] for a slightly different nonimaging 
concentrating compound parabolic trough solar collector 
operating under slightly different atmospheric conditions. 
It is interesting to note that the larger capacity solar 
collector had the best instantaneous efficiency and also 
the highest life-cycle fuel savings. In all three cases 
the optimum incident acceptance angle was found to be 
18.04 - 18.05 degrees and the optimum truncation angle 


was found to lie within the range of 184.55 to 190.0 degrees. 
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Due to the fact that the objective function is weakly 
linked to the design capacity of the solar collector, the 
final daily heat gain 1s somewhat higher than the minimum 
set. If a stronger link were to be established between 
the solar collector capacity and the objective function, 
the collector design would be driven closer to the stated 
Minimum daily average heat gain. Also, if the convergence 
criteria is tightened, the optimizer will take longer to 
reach the minimum but will reach a final design where the 


average daily heat gain is closer to the minimum set. 
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APPENDIX A 
PESORel os | Se ROGRAM RESULTS 
This appendix contains the results of the four test 

programs that were used to develop the DESOP numerical 
Optimization program. For each design the initial design, 
the true optimum design, the DESOP results using default 
control parameters, and the DESOP results using adjusted 
control parameters are given. The four test programs may 


be found in Appendix F. 
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DESOP 


Pent. Suboroqram BANNA 


Initial Design: 
Design Variables: 


wal} = —1.2 
X(2) = 1.0 
Segeceive Function: 
Obj = 10.8 
Oseg = 10.8 
Side Constraints Violated: 
N/A 
Constraints Violated: 
N/A 


DESOP Results: 


Default Control Parameters: 
Design Variables: 


X(1) = Q.768 
ye) = 0.578 
Objective Function: 
Goji = 4.05 
Opj = 4.05 
Seae Constraints Violated: 
N/A 
Constraints Violated: 
N/A 
# of Function Evaluations: 
96 


NOTE: 


constraints. 


The BANNA subprogram has 


iio oi Ci nA 


True Optimum: 
Design Variables: 


X(1) = 1.00 
<i = 1200 
Objective Function: 
Obj = 4.00 
Opj = 4.00 
Side Constraints Violated: 
N/A 
Constraints Violated: 
N/A 


Adjusted Control Parameters 
Design Variables: 


KCL) = 0.791 
xX{(2) = 0.606 
Objectuve Funceion: 
Obj = 4.047 
Opj = 4.047 
Side Constraints Violated: 
N/A 
Constearnes Violatecd: 
N/A 
# of Function Evaluations: 
91 


2 design variables and no 





| 
| 
| 


DESOP TEST PROGRAM 


ANALIZ Subprogram : 


Initial Design: 
Design Variables: 


Rie) = 1 
X(2) = 1 
X¥(3) = 1 
X(4) = 1 
Objective Function: 
Obj = 3l 
Opj = 31 
Side Constraints Violated: 
N/A 
Constraints Violated: 
None 


DESOP Results: 


Default Control Parameters 
Design Variables: 


X(1) = 4.72E-02 

X(2) = 0.998 

X(3) = 1.98 

X(4) = -1.00 
Objective Function: 

Obj = 6.088 

Spy = 6.093 
Side Constraints Violated: 

N/A 


Constraints Violated: 
G(3) = 0.000527 

# of Function Evaluations: 
392 


Resen-Suzuki 


True Optimum: 
Design aes 


X(1) = 
M2) = 0 
X(3) = 2.0 
X(4) = -1.0 
Objective Function: 
On je = 100 
Opj = 6.00 
Side Constraints Violated: 
N/A 
Constraints Violated 
None 


Adjusted Control Parameters 
Design Variables: 


Mt) = —-5-167E-03 
X(2) = 1.019 

KiC3) = 1.999 

X(4) = -0.9951 


Objective Function: 
Obj = 5.9998 


Opae= 6.007 
Side Constraints Violated: 
N/A 


Constraints Violated:. 
G(3)= 0.00232 
# of Function Evaluations: 
306 


NOTE: The Rosen-Suzuki subprogram has four design variables 


and three constraints. 
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Paes, Subcroaram : 


PeSsOr wlfot © ROGRAM 


T/VAR 


Paacial Design: 

Design Variables: 

X (1) il 
X(2) = 
X(3) = 
X(4) = 
reCS) i= 
X(6) = 
C 


I 


m{ } 
Side Cons 

None 
Segective Function: 

Obj = 7-203 

Opj = 1007 
Constraints Violated: 

G(l) = Ll 


= 
1 
it 
1 
iL 
a 
a 


raints Violated: 


DESOP Results: 


Default Control Parameters: 
Design Variables 


mil) = 1,42 
mee2) = 0.515 
m3) = 0.376 
X(4) = 0.844 
mes) = 2.13E-02 
X(6) = 9.84E-03 
m7) = 0.638 
Side Constraints Violated: 
None 
Objective Function: 
Obj = -142 
Opj = -142 
Constraints Violated: 
None 
# Of Function Evaluations: 
302 


(D 


Meews:, ine T/VAR subprogram nas s 


Genser aint. 


True Optimum: 


Design Variables: 
X (1) 
m2) 
3) 
xX (4) 
X(5) 
X(6) 
X(7) = 

Side Constraints Violated: 
None 

Objective Function 
©by = ~190 
Opj = -190 

Constraints Violated: 

None 


| 
SOOfrP COW 


Adjusted Control Parameters: 


Design Variables: 


X(1) = 2.69 
X(2) = -7.60E-04 
Moo) os 1. 04e—0 2 
M2 je =  bac2 
X(5) = 2.89E-03 
X(6) = 2.91E-03 
X77) = 0.213 
Side Constraints Violated: 
1 
O©byjective Function: 
Obj = -178.6 
Gog. L719. 1 
Constraints Violated: 
None 
TeOeee 1ceELON Evaluations: 
808 
Vemmeesl ene variables ana one 
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APPENDIX B 


NISCO DESIGN RESULTS 


| DESIGN SPECIFICATIONS 


A. 


GEOGRAPHIC 


Location 
Solar position and 
intensity levels. 


Wall azimuth angle 
Cloud Cover 

Vapor Pressure 
Tamb 


COLLECTOR MATERIALS 


RECEIVER 
Absoptivity 
Thermal Emissivity 
Reflectivity 


REFLECTOR 
Reflectivity 


COVER 
Specular Absorptance 


40 degrees North Lattitude 

ASHRAE Handbook of Fundamentals 
standard solar radiation design, 

Ref. 

0 deg. 

0 % 

3 mm Hg 

ASHRAE Handbook of Fundamentals 
Gaily norms for San Francisco, Ca. 
Ref. 


Oxidized Copper - 
0.93 
0.40 
0.07 


Vacuum deposited Aluminum on resin 
0.89 


Double strength window glass 
ASHRAE standards for D.S. window 


Specular Transmisivity glass. Ref. 


Absorptance av. 
Transmittance av. 
Reflectivity av. 
Thermal Emissivity 


COOLANT | : 


Specific gravity 
Specific heatr 
ViscositytT 

Inlet Temperature 


0203 
0.80 
Ged 
0.94 


Therminol 55 


Oe o7 

4.9E~4*Tr+0.4036 Btu/(lbm F) 

bree goa 4*(=0.0534*7r+32.3) lbm/ (ft sec) 
100 deg. F 


tT Values given as a function of average receiver temperature, 


Tr, in degrees F. 
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ECONOMIC 


Life-Cycle 

Fuel Cost Ree Gas) 
Annual Fuel Inflation Rate 
Monetary Inflation Rate 
Gollector Cost 


20 years 
8.14 $/10° 
Oil: 

0.10 

22 So Et: 


64 


Btu 


Ref. 


d 





_ NiSCcO BESIGN = 10000 Btu/DAY SOLAR COLLECTOR 


INITIAL DESIGN 
Design Variables: 


—— ee 


65 


Incident acceptance half angle 25.00 deg 
Truncation angle 230.0 deg 
Receiver radius 0.41 in 
Collector length 20-00 -£t 
Coolant mass flow rate 108.7. 2bm7 ir 

Design Features: 
Concentration ratio 2.20 
Gollector aberature area 9.56 ft? 
Collector depth 4. 2m 
Coolant velocity Os taert/sec 
Average daily heat gain 11,400 Btu 
Maximum coolant temperature 1357 deg F 
Instantaneous collector efficiency 0270 
Initial cost $203.00 
Life-cycle cost savings $7.78 VAS? Btu 

FINAL DESIGN 

Design Variables: 
Incident acceptance half angle 18.05 deg 
Truncation angle 190.0 deg 
Receiver radius Oe 52)FIn 
Collector length to oe 
Coolant mass flow rate 1243 Salam ae 

Design Features: 
Concentration ratio LoS 
Collector aperature area 8.48 ft* 
Collector depth Oe en 
Coolant velocity Opi let e7sec 
Average daily’ heat gain 10,600 Btu 
Maximum coolant temperature 2s deg F 
Instantaneous collector efficiency 0.74 
Initial cost $159.00 
Ieere-cyecle cost Savings $8.15 /10° Btu 





PeENITIAL 


DESIGN: 


Design variables: 


Incident acceptance half angle 
Truncation angle 

Receiver radius 

Collector length 

Coolant mass flow rate 


Design features: 


Soncentration ratio 

Collector aperature area 

Collector deoth 

Max. coolant temperature 
Instantaneous collector efficiency 
Average daily heat gain 

mattlal cost 

Life cycle savings 


PINAL DESIGN: 


Design variables: 


Incident acceptance half angle 
Truncation angle 

Receiver radius 

Collector length 

Coolant mass flow rate 


Design features: 


Run time - approx. 


Concentration ratio 
Collector aperature area 


iSellector depth 


Coolant velocity 

Max. coolant temperature 
Instantaneous collector efficiency 
Average daily heat gain 

Paitial cost 

Life-cycle savings 


5S nours 
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NISCO DESIGN - 30000 Btu/DAY SOLAR COLLECTOR 


Po 
230.0 

1.00 
20.0 
100.0 


eo 
Peles 

M2 Ore 
366.59 
0.64 
83,400 
$1,628.90 
$7.54 


18.04 
64.25 
£2.28 
26.43 
£25523 


1.54 
27.395 


6.04E-02 
i294 
Or 
34,473 


9909.98 
98.19 


ft: 
lbm/hr 


ft2 
aati 


Btu 


7 02 Btu 
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APPENDIX’ € 


DESOP USER'S MANUAL 


Bes INTRODUCTION 

The numerical optimization package of programs, currently 
consisting of OPCON and DESOP, provide the capability for 
finding the optimum design of a system mathematically modeled 
using multiple variables, on the Hewlett-Packard 9845A desk- 
top computer. OPCON iS an executive program which provides 
the user with the following: a primary point of contact with 
the computer from which to access the optimization program 
DESOP, a standard formatted input for design variables, side 
constraints on the design variables, optimizer control varia- 
bles, and organization of the optimization process. DESOP 
1s a numerical optimization program, which when coupled to 
a user supplied design analysis program, will optimize the 
design. DESOP will allow the user to monitor the optimiza- 
tion process as it is taking place. After monitoring the 
optimization process, the user may chooseto change the 
optimizer control variables and/or his design starting 
point to more efficiently or more accurately reach the 
design optimum. After an optimization has been performed, 
DESOP will reload the OPCON program and return the user to 


the OPCON program. The OPCON program will then offer the 
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user the following three choices: to optimize the same 
design, optimize a different design or terminate the 
program. 

Figure ll illustrates the overlay method used in 


loading the OPCON and DESOP programs. 
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PROGRAM OVERLAY STRUCTURE MASS STORAGE 
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Figure ll. Basic Program Relationships 
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Peet, CONCEPT OF NUMERICAL OPTIMIZATION 

Consider the problem of designing the cantilevered beam 
shown in Figure 1. The design task may be broken down into 
three major oarts. First, the objective of the design must 
be determined, which in this case is to minimize the weight 
of the beam required to support the concentrated tip load P. 
Second, any physical constraints that may effect the design 
fms. be determined. Thirdly, any limits which exist on the 
design variables must be determined. The design problem 
may then be reduced to a system of equations as follows: 

Minimize the volume (V) 

V = B-H-L 
Suoject to: 


Bending stress (c,) 


Ga 
Oo, = amar 2 se20000 pst 


Shear stress (v) 


eo 2 ee 
Moma 10000 psi 


Mer leece7On wumeer load (4) 
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P=10000 





Figure 1. Cantilevered Beam Design Problem 
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With geometric constraints such that: 


Ceo eens 0.5 
ie Ome it 20.0 
babe OO 


At this time the following definitions are introduced: 

Gbjective Function - The: parameter which is to be 
minimized or maximized during optimization. The objective 
function always occurs on the left side of the equation 
unless it 1s also used as a deSign variable. An objective 
function may be either linear or nonlinear, implicit or 
explicit, but must be a function of the design variables. 

Design Variable - Any parameter which the optimiza- 
tion process 1s allowed to change in order to improve the 
design. Design vals appear only on the right hand 
Side of equations in the analysis program. 

Inequality Constraint - Any parameter which must 
not exceed specified bounds for the design to be acceptable. 
Constraint functions always appear on the left side of 
equations. A constraint may be linear, nonlinear, implicit 
or explicit, but must be a function of one or more design 
variables. 

Equality Constraint - Any parameter which must equal 
a specified value for the design to be acceptable. The 
Same tules apply to equality constraints as inequality 


Sonstraints. 
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Side Constraint - Any upper or lower bound placed 
upon a deSign variable. Side constraints are usually not 
included in the system of equations that comprise an 
analysis program. Instead they are usually included as 
Peiceeor the data input to the optimization program. 

Analysis Code - The system of equations utilizing 
the design variables which are used to calculate the objective 
mmetiOn and the constraints of a particular design problem. 

The general optimization problem may thus be stated 
mathematically as: 
Find the set of design variables X. Where: 1 = Wie). ee 
Mauch will: 
Minimize the objective function (0b}3) 
Obj = £(X) 
Sooyect Eo: 
Inequality constraints (G) 


G, (es 30 eee, ia 


Equality constraints (H) 


H, (X) = @ oe ele ae) 
Side constraints 

xt < x, < xl tee oe Th 

hs a: as 


Returning to the cantilevered beam problem, it may be 
Seaeea in the standard format as follows: 
Memext!) = 3, <(2) = H, and Obj = Vol = B-H-L 


Then Minimize Cbj = Vol 
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SUDZeECE tO: 


oe 
ee ee. < 
Geil) 50000 LO 
eee 
G(2) = 999 ~ 1 = 
eee 0 
G(4) =2-10< 0 


With side constraints: 
il 


AG =O 265 
ce 50 
x(2)* = 1.0 
a2) c  =200 


Zt is thus fairly simple and straightforward to perform 
an analysis on a particular beam for a given B and H. 
Successive analyses may be performed on the cantilevered 
beam by solving the above system of equations. It is de- 
Sirable to automate the succesSive solutions and to direct 
the solutions such that each solution 1s a better design 
than the last. One approach for doing so, and the one 
used by DESOP is to proceed as follows: Start with initial 
values for B and H. Solve the above set of equations to 
Bema the objective function Obj and to see if any constraints 
are violated. <A pseudo objective function is created to 
represent designs when constraints are violated. Ifa 


a 


Bemoecralnt is viclatec, a penaltv is added to Obj to form 


1 


pepeemalized Objective function Onj. The gradient of the 
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penalized objective function at the initial design may be 
found by taking the first partial derivative of Opj with 
respect to the design variables. The gradient of the 
penalized objective auncsdee G@etines thesdirectiensot 
peeepest ascent. In the case of the cantilevered beam, it 
is desired to minimize the objective function; therefore, 
the greatest improvement in design may be achieved by moving 
in the negative gradient, or steepest descent direction. 
From the initial design point a search is performed in the 
steepest descent direction for the minimum value of Opj in 
that direction. At the new minimum, the gradient of the 
penalized objective function 1S again determined and a 
search 1s performed in a conjugate direction until a second 
minimum is found. Successive iterations are performed 
until the gradient is found to be zero or each successive 
iteration produces a sufficiently small change in Opj such 
that for all practical purposes the minimum has been found. 
At this time the penalty function is increased. If the 
design is in a region where there are no constraints 
violated an increase in the penalty function will not 
Smange the value of Opj. I£ on the other hand the design 
2S in an infeasible region where there are one or more 
Seastraints violated, Crj will be increased, and the 

search for a new minimum will commence. If the minimum 


Seemtne Obtective function exists in the infeasible region, 
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the minimum value for the objective function in the feasible 
region will be approached from the infeasible region as the 
penalty function 1s increased. The design improvement 
process will terminate when a zero gradient is found or 
Succes3ive iterations produce a sufficiently small change 
in the value of Opj and an increase in the penalty function 
causes no change in Op}. 

The minimum thus reached by the optimization process is 
a Minimum with respect to the penalized objective surface 
immediately surrounding the final design point. The opti- 
mization process cannot distinguish between local and global 
minimum points. It is thus good engineering practice to 
run several optimizations for a particular design problem 
Peom several different initial design points. If optimiza- 
tions performed from different initial design points converge 
on the same minimum point, that point 1s probably a global 
minimum. If on the other hand two or more minimums are 
found, there may be local minimums located in the design 
Space being considered and care must be taken to find the 


global minimum. 


BF Gime DESOP NUMERICAL OPTIMIZATION PROGRAM 

The DEsktop Sequential unconstrained minimization 
Meemmicgue Optimization Program was developed using the 
Meee ODtIimization apvroach cutlined in Section II.A. 


= 


Meer tie O2OGram 1S included in Appendix E. The major 


~ 
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program structure is shown in Figure 2. The following 
discussion will refer to Figure 2 and describe the major 
features of the program. 

7. Basic Program Execution 

| The DESOP program begins execution when it is 

loaded, linked to the user's analysis subprogram, and the 
program is instructed to run by the OPCON program. The 
above actions are automatically performed by the OPCON 
program. DESOP 1s loaded into the computer by an overlay 
process. Therefore no variables can be directly transferred 
between the DESOP and OPCON programs. DESOP begins execu- 
tion by reading the optimizer control variables and the 
design variables that were input using the OPCON program 
and saved to a mass storage device. The program then sets 
Icalec equal to one and evaluates the objective function 
and constraints at the initial deSign point. Icalce is a 
flag provided the user to key user specified output on the 
initial and final design analysis. DESOP will provide the 
user with a hard copy output of the initial design variables, 
the value of the constraints, the objective function and 
the penalized objective function. The user then has the 
@@eron Of continuing with the DESOP program to optimize 
his analysis subprogram or to return to the OPCON program 


GO change one or more of the input parameters. 
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INPUT X5 
| EVALUATE Obj & G 


FIND THE GRADIENTS 
OF Obj AT X 


CALCULATE THE 
FLETCHER REEVES | 
SEARCH DIRECTION 






















CALCULATE THE 
Srebrhot DESCENT 
SLARCH SOLTRECTION 


ESTIMATE a FOR MIN. Obj 
IN 1-D SEARCH 











FIND a FOR MIN. Ob} 

BY EXPANDING SEARCH, 
METHOD OF GOLDEN SECTION, 
& POLYNOMIAL INTERPOLATION 


ee 


MOVE TO NEW X 
CHECK CONVERGENCE OF AOQbj 


INCREASE PENALTY FUNCTION 
RECHECK CONVERGENCE 


Paguce 2. DESOP Flow Diagram 





Proceeding with the optimization, there are two 
major loops in the optimization program. The outer loop 
M@emeases Ehe penalty function when the inner loop's con- 
vergence criteria have been met. A convergence test is 
then performed by the outer loop. If the convergence 
criteria 1S met, the optimization process is considered 
finished. If the convergence criteria for the outer loop 
1s not met, program execution 1s returned to the inner 
loop. The inner loop performs successive iterations 
Searching for the minimum of the penalized objective 
function with no increase in the penalty function taking 
place. When the inner loop's convergence criteria have 
been met program execution is transferred to the outer loop. 

Execution of the program while in the inner loop 
proceeds as follows: First, the gradient of the penalized 
objective function is calculated by subroutine GRAD. The 
program then computes a search direction using either a 
steepest descent method or the method of conjugate directions 
developed by Fletcher and Reeves [Ref. 1]. Once a search 
direction is established the optimizer attempts to locate 
the minimum value of the penalized objective function in 
m@@ensearch direction. This process if referred to as the 
One-dimensional search and is illustrated in Figure 3. 

Mne efficiency and accuracy to which the one-dimensional 


Search for the minimum cf the penalized objective function 
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is accomplished, 1s the key to successful sequential un- 
constrained minimization technique numerical optimization. 
The one-dimensional problem may be expressed in terms of 
the penalized objective function, Opj, and the amount of 
movement, «, in the search direction, S. First the slope 
of Opj with respect to + is calculated. An initial "guess" 
of how far to move 1S made uSing subroutine ALPGES. The a 
which corresponds to the minimum value of Opj in the one- .- 
dimensional search 15 then calculated using subroutine 
ALPBND and subroutine QUEBIC. The minimum value of Opj 
thus found is then compared to the previous value of Op} 
for convergence uSing subroutine CONVRG. If convergence 
mmo met, execution returns to the start of ‘the inner 
loop. If convergence 1S met, execution returns to the 
Outer loop. 

When the convergence criteria have been met for 
Seieaeeie inner anc outer loops, the program proceeds to 
Set ICALC to three as a flag for user generated output for 
the final design. DESOP then provides the user with a 
hard copy output of the final design variables, objective 
function value, penalized objective value, constraint 
values, the number of inner loop iterations, the number of 
memes the analvsis suborogram was called and the final value 
Meee Penalty function. “fhe OPCON program is then overlaved 
Over the CESOP program and program execution is returned to 


eae OPCON program. 
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Pewee binding the Search Direction 
(ieee p In ftinadumg the search darection, S, 
Memeo determine the slope of Opj at the present design 


Peeme. The forward Finite difference method is used where: 


e + * i a * 

ee F(X; + AX;) - F(X;) © * 

Ax ok. ia i 
aL J 


2S ae eeLON. 
As dF / Gx, gives the direction of positive slope, the search 
direction 1s the negative of JES OX, . The first search 1s 
perfcrmed using the steepest descent as found above using 


the following relation: 


ee = xX. + as. 
le aL au 


where alpha is the distance moved in the S$ ike oe When 
a Minimum is obtained along the direction of steepest 
descent, a new Fletcher-Reeves conjugate search direction 
[Ref. 1] is calculated at the new Design point using the 


following relations: 








Se as. 
ox it 
Ndv 
v oa 
‘sal aX. 
Be i=l a 
Ndv JF ,- 
y { xT 
i=] 


where the prime denotes values for the present iteration 


and the non-prime variables indicate values for the previous 
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iteration. A one-dimensional search is then performed in 
the new search direction. Searches are continued using 

the conjugate direction method for Ndv + 1 iterations, 

where Ndv is the number of design variables. The search 
process 1s then restarted using the steepest descent method. 
The reason for incorporating the conjugate search method 

is that the steepest descent method when traversing a design 
surface with a curved valley will tend to zigzag from one 
Side of the design surface valley to the other making very 
little progress as 1S illustrated in Figure 4. The conju- 
gate direction method is much more efficient in traversing 
such a design surface. However, as the conjugate direction 
method is additive upon previous searches, it has a tendency 
to decrease in effectiveness with each successive search 


owing to the accumulation of numerical "noise." For that 
reason the search process 1S restarted with the steepest 
descent method every Ndv + 1 iterations, or when the con- 
jugate direction predicts a positive slope. The search 
direction is normalized to avoid inaccuracies caused by 
Gemerical ili-conditioning. 
Pee, Estimating an [Initial Value for Alpha 

The initial estimate for alpha is made in the 
following manner: First, the slope of Opj in the search 
erection is calculated as the sum of each of the products 
Seeene Gradients times tne search direction. Then the 


meee Of Op] in the search direction is divided by tne value 
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of Opj. This value is then multiplied by an improvement 
percentage in Opj. This first estimate is then applied to 
a series of conditional tests to determine the'validity of 
the estimate with respect to the slope of Opj and the magni- 
tude of the design variables. Lastly, the estimate for 
alpha is checkec to see if it violates any side constraints. 
If it does, the value of the estimate for alpha 1s reduced 
until the side constraints are no longer violated. 

4. Calculating Alpha 

The calculation of alpha is the most critical al- 

Sorachm in the DESOP program in providing reliable optimizer 
Meeration. The ability to accurately and efficiently find 
the minimum of the penalized objective function in the one- 
dimensional search affects directly the operation of the 
optimizer. Figure 5 illustrates the zigzag phenomenon 
which occurs if alpha is not accurately found. The zigzag 
phenomenon is caused by the fact that the optimizer in per- 
forming the forward finite difference for calculating the 
search direction perturbs the design vector a very small 
amount. As such the optimizer can only "see" the design 
surface that is immediately adjacent to the design point. 
Therefore, if the minimum is not found in the one-dimensional 
Search, the optimizer will converge very slowly on the 
minimum. 

There are tow major secticns to the ALPBND subroutine. 


Mae Lirst section attemots to find the minimum value of Opj 
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Figure 5. The Zigzag Phenomenon 
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uSing an expanding search technique. The first move is 
the amount predicted by the ALPGES subroutine. If the 
minimum is not bracketed by the first move, additional 
moves are made. Each move is larger than the last. The 
Size of the move 1S increased each time by an amount equal 
to the size of the last move divided by the lower Golden 


Peetion fraction, where the Golden Section fractions are: 


3 - 75 


The lower Golden Section fraction, F is used so that the 


L’ 
interval will be consistent with the Golden Section search 
in the second section of the ALPBND subroutine. The expand- 
ing search is continued until the minimum value of the ob- 
jective function has been bracketed. 

Once the minimum is bracketed, a Golden Section 
search is performed to reduce the bracketing interval on 
the minimum by an amount such that when the two end points 
of the interval are taken with two points internal to the 
interval and a cubic is passed through the four points, the 
Cubic will accurately predict the minimum of the penalized 
objective function. Himmelblau in [Ref. 2] states that the 
Golden Section search method of reducing the interval 
eround the minimum of Op} is the most effective of the 


reducing techniques studied. Golden Section search is based 
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on the splitting of a line into two segments known in 
ancient times as the "Golden Section." The ratio of the 
whole line to the larger segment is the same as the ratio 
of the larger segment to the smaller segment. The two 
Golden”Section fractions are emploved to split the interval 
bracketing the minimum as shown in Figure 3. Once the 
interval has been split, the two values of Opj corresponding 
to the internal points are compared to find the larger of 
the two. The internal point with the larger value of Opj 
will become the new end point for the interval, the remain- 
ag Meer ior point will by the fact that it was determined 
by a Golden Section fraction, be equal to the point deter- 
jmmea by the other Golden Section fraction. Thus, only one 
new point must be calculated to continue the Golden Section 
search. The search is continued inch manner until the 
vertical separation of the two end points with respect to 
the interior points is less than one percent. The four 
values of the penalized objective function corresponding to 
the four Golden Section search points are then sent to a 
Be@erc interpolator. The cubic interpolator will return a 
value for alpha to predict the minimum of the penalized 
objective function, and the nin ie Sueene Clore. funceton 
that the interpolator has created. The subroutine ALPBND 
will then test the predicted minimum with the minimum found 


at the predicted alpha. If there is less than a tenth of 
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one percent difference between the two values of the objec- 
tive function, the point predicted by the cubic interpolator 
will be accepted as the minimum and program execution will 
meeurm to the main program. If the predicted minimum is not 
sufficiently close to the minimum at the predicted alpha, 
another Golden Section search will be performed to reduce 
the interval and better localize the minimum. The four 
points from the reduced interval will then be sent to the 
cubic interpolation subroutine. This process will continue 
until either the test for the minimum is positive or the 
interval has been reduced to less than 1LE-12. Program exe- 
cution will then return to the main program. 
PemeoubDroutine QUEBIC 

Subroutine QUEBIC 1s used to estimate the alpha at 
which Opj iS a minimum based on four point cubic interpola- 
tion. If the function more closely resembles a quadratic 
than a cubic, a three point quadratic interpolation 1s 
performed using the three points which bracket the minimum. 
If the predicted minimum 1s outside the interval spanned 
by the two end points again a quadratic interpolation is 
performed. If the minimum still lies outside the two end 
points, the analysis returns to subroutine ALPBND, the 
inverval bracketing the minimum 1s reduced, and program 


Seeecution returns to QUEBIC. 
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Pe oOmve agence Or che Penalized Objective Function 


The penalized objective function is tested for con- 
vergence at the end of each inner loop and again at the end 
of the outer loop in the main program. Convergence is 
tested by calling subroutine CONVRG. There are two criteria 
used for testing for convergence. The first tests the rela- 
tive difference of the value of Opj from the present iteration 
with the value of Opj from the last iteration. The second 
method tests the absolute difference of the two values. The 
second method is emploved for cases when the value of the 
penalized objective function approaches zero. When conver- 
gence has been met on two successive iterations, the penalty 
Function is increased by an amount specified by the user in 
the executive OPCON program. The penalized objective function 
1S again tested for convergence. If convergence 1s still 
met, the optimizer considers the present value of the pena- 
lized objective function to be a minimum, noting again that 
numerical optimization programs cannot differentiate between 
local and global minimums. 

fee Ene Penalty Function 

The purpose of the penalty function is to increase 
the value of the objective function when the design is in 
an infeasible region. The infeasible region is that region 
where one or more Gesign constraints are violated. When a 


Bemsctralint iS violated, the value of the particular constraint, 





Gyr 1S poSitive. The objective function is then penalized 


as follows: 


Ong — Ob) Ga 
where: 

Ro> a murtiplication conmstant 

Peoean CxPOnentc Constant 
This type of penalty function, one where the penalty 1s 
applied aiter the design leaves the feasible region, 1s 
known aS an exterior penalty function. The exterior type 
of penalty function was chosen over other types, such as 
the interior or extended interior penalty function. Ifa 
function is discontinuous within the design space being 
studied, numerical difficulties may be encountered which 


make performing an optimization of the design difficult. 
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evo OF THE DESOP NUMERICAL OPTIMIZATION PROGRAM 

Once the analysis subprogram has been written and SAVED, 
the user may perform a numerical optimization of the design 
analysis. Figure ll shows the basic relationships between 
the OPCON and DESOP programs and the user supplied ANALIZ 
Subprograms. Note that the DESOP program overlays the OPCON 
program after all data has been input. To begin, LOAD OPCON 
and press run. Through a series of self-explanatory menus, 
the user will be prompted to input control variables, design 
variables and execute the DESOP program. The following 
menus appear in the OPCON program. 

Ll. INTRODUCTION - A brief introduction stating the 
purpose of the OPCON program. 

2. NEW OR EXISTING ANALYSIS PROGRAM - If this is the 
first time that a particular analysis subprogram has been 
run on the tape or disk being used, or if the files from 
a subsequent run have been deleted, the user must enter 
the response for a new analysis subprogram. For existing 
programs, OPCON will read data from the existing data 
file. The data read will be that from the subsequent run 
of the particular analysis subprogram. 

3. NAME THE ANALIZ SUBPROGRAM - The user will input 
the name under which the analysis subprogram he wishes to 
use was saved. If this is a new analysis subprogram 
OPCON will create files for saving the data input during 


the execution of OPCON. 
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4. OPTIMIZER SELECTION - At the present time there is 
only one optimization program available. It is hoped that 
at some future date additional numerical optimization 
programs will be added to the optimization package. OPCON 
has been developed to interface with multiple numerical 
optimization programs. Details for linking other numerical 
optimization programs to OPCON are included as comments in 
the OPCON listing which may be found in Appendix D. 

5. INPUT NDV AND NCON - For a new analysis subprogram 
the user will be asked to input the number of design varia- 
bles, NDV, and the number of constraints in the analysis 
Subprogram, NCON. 

Si. INPUT CHECK OF COMMON CONTROL VARIABLES ~- A menu of 
the optimizer control parameters common to DESOP, Feasible 
Direction type and other future optimization programs is 
displayed. The menu displays the variable name, its 
minimum, maximum, default and present value. The variables 
displayed are: 

NDV - The number of independent design variables 
used in the analysis code. 

NCON - The number of constraints in the design 
analysis subprogram. 

IPRINT - A print control used in the optimization 
program to display intermediate results. Positive values 


entered will print on the CRT. Negative values entered 
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will print on the thermal printer. If a zero is entered 
there will not be a hard copy output of the initial and 
final results of the optimization program. IPRINT = 0 is 
to be used when debugging an analysis subprogram to 
conserve thermal paper. 

i een nee 1) creer ao p . Lieaat LON 
information. 

gee olin abeve pilus the objective 
function and penalized objective function on each iteration. 

foes Cerne above plus=cthe Constraine 
values, search direction vector and move parameter alpha on 
each iteration. 

Poa rte above Pius Gradient informacion 
on each iteration. 

jo soni noes Mis cach Proposed, design 
vector and the penalized objective function during the one- 
dimensional search on each iteration. 

Tom, Debugging ‘alcmror Cptimizer devellop= 
ment. DESOP will pause after each major operation is 
performed during the optimization process. 

DELFUN - The minimum absolute change in the objec- 
tive function to indicate convergence of the optimization 
process. 

DABFUN - The minimum absolute change in the objec- 
tive function to indicate convergence of the optimization 


PlOcCess. 


oD 





ITMAX - The maximum number of inner loop (uncon- 
strained minimizations) without increasing the penalty 
Punction. 

ICNDIR - The conjugate direction restart parameter. 
Every ICNDIR inner loop iterations a steepest descent search 
1s performed. It is recommended that ICNDIR be set equal 
mOeNDV + 1. 

FDCH - The relative change in the design variables 
for calculating finite difference gradients. 

FDCHM - The minimum absolute step in finite difference 
gradient calculations. 

ABOBJ1 - The expected fractional change in the ob- 
Neeerve function for the first estimates of the step size 
to be taken in the one-dimensional search. 

ALPHAX - The maximum fractional change in any design 
variable for the first estimate of the step size to be taken 
in the one-dimensional search. 

ye. PNeEUD CHECK OF SfHESDESOP CONTROL PARAMETERS - A menu 
of the optimizer control parameters for the DESOP optimiza- 
tion program is displayed. The menu will display the variable 
name, its maximum, minimum, default and present value. The 
variables displayed are: 

IRMAX - The maximum number of times that the penalty 
parameter will be increased. 


RZ - The starting value of the penalty parameter. 
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execution beginning at line OPT1. The storing and reloading 
process allows the OP program to be called and run as a 
separate program for debugging purposes. At the completion 
of the optimization process, the optimizer program will 
reload OPCON and return execution to the OPCON program. A 
menu will be displayed welcoming the user back to the OPCON 
program giving him the following options: to optimize 

again using the same ANALIZ subprogram, to optimize again 
uSing a different ANALIZ subprogram, or to terminate the 


program. 
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APPENDIX E 


DESOP Program Listing 


In the development of the DESOP program, Refs. 2, 7, 8, 


9, and 10 were used. 
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